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ABSTRACT 

m 
o , 

O , We describe DEVA, a multistep AP3M-like-SPH code particularly designed to study 

galaxy formation and evolution in connection with the global cosmological model. This 
code uses a formulation of SPH equations which ensures both energy and entropy con- 
servation by including the so-called V/i terms. Particular attention has also been paid 
to angular momentum conservation and to the accuracy of our code. We find that, in 
order to avoid unphysical solutions, our code requires that cooling processes must be 
^ . implemented in a non-multistep way. 

?N ; We detail various cosmological simulations which have been performed to test our 

^ ', code and also to study the influence of the S/h terms. Our results indicate that such cor- 

I rection terms have a non-negligible effect on some cosmological simulations, especially 

I on high density regions associated either to shock fronts or central cores of collapsed 

' objects. Moreover, they suggest that codes paying a particular attention to the im- 

■ plementation of conservation laws of physics at the scales of interest, can attain good 

accuracy levels in conservation laws with limited computational resources. 

o 

H 

I Subject headings: galaxies: formation - galaxies: discs - hydrodynamics - methods: 

^ ' numerical. 
> 

^- 

^ . 1. Introduction 

In the last few years most cosmological parameters have been determined up to a few percent. 
The values of Otot, ^A, ^^baryon and Hq can now be constrained with an unprecedent degree of 
accuracy (see, for example, Lahav 2002; Lahav et al. 2002; Netterfield et al. 2002; Spergel et 
al. 2003, and references therein). The next challenge to cosmologists is to test the predictions of 
cosmological models at a few hundred kpc scales. It turns out that these are just the relevant 
scales involved in galaxy formation and evolution. Galaxy formation and evolution are intriguing 
open questions whose resolution in connection with the global cosmological model will very likely 
advance considerably in this decade. Even though the field is at its beginnings, the use of numer- 
ical methods to study how galaxies are assembled within a cosmological scenario from the field of 
primordial fluctuations, seems a convenient approach. The main advantage of these approaches 
(i.e., self-consistent hydro dynamical simulations), is that the physics is introduced at a very gen- 
eral level, and the system evolves as a consequence. We can follow the evolution of the dynamical 
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and hydrodynamical properties of matter in the Universe; galaxy-like objects (GLOs) appear as 

a consequence of this evolution. And, so, the building-up of objects (cosmic network structure 
formation at high z, collapse, interactions, mergers, accretions), as well as their hydrodynamical 
consequences (instabilities, gas infall from halos to discs at hundred of kpc scales, gas inflow along 
discs at tens of kpc scales, shocks, cooling, piling- up of gas necessary for star formation), can be 
accurately followed. We get not only the properties of objects at any z, but also an insight into 
the physical processes responsible for their formation and evolution. Moreover, numerical hydro- 
dynamical simulations using particles permit very convenient comparisons of GLOs that form in 
simulations with observational data. Simulations directly provide us, at each z, with the structural 
and dynamical properties of each individual GLO (position and velocity of each of its particles, 
gas density and temperature of each of its baryonic constituents) and with their individual star 
formation rates histories (SFRHs). Chemical abundance and spectrophotometrical data are the 
current standard to compare models of galaxy formation. It is expected that the next generation 
of astronomical facilities will make possible a new science: mass measurements for distant galaxies 
(sec, for example Vcrhcijcn, Bcrshady & Andersen 2002). GLOs formed in numerical simulations 
arc particularly suited to be compared to this new kind of data. 

Pre-preparcd numerical experiments are adequate to describe in detail a particular phase of 
the formation or evolution of galaxies (for example, merger events or orbital motions of satellites 
within halos), from initial conditions set by the experimenter. These initial conditions try to model 
conditions that would have arisen along the evolution of the systems under consideration. They 
are useful to study basic aspects of the physical processes relevant to evolution. For example, 
the works by Barnes, Hernquist and Mihos (Barnes 1989; Barnes & Hernquist 1991; Barnes 1992; 
Barnes & Hernquist 1992; Mihos & Hernquist 1994, 1996), which have fundamental importance 
to understand the role played by mergers in galaxy evolution, have been carried out with this 
technique. However, in pre-prepared simulations, contrary to the self-consistent approach, the 
process under consideration is studied in isolation, and not in connection with the other relevant 
processes involved in galaxy formation and evolution (already mentioned) that, moreover, could 
interact among themselves in a non-trivial way. 

These considerations stress the ability of self-consistent hydrodynamical simulations as a tool 
to learn how galaxies form from the field of primordial fluctuations and evolve into the objects 
we observe today. To properly handle this problem, a numerical code has to allow for enough 
mass, time and space resolution as well as a convenient dynamical range. They should be as fast as 
possible and with memory requirements within the current computer capabilities. A very important 
issue when designing a numerical code to study galaxy formation and evolution, is to make sure 
that conservation laws are accurately verified, and, particularly, i), that angular momentum is 
conserved at the scales relevant to disc formation; otherwise, galaxy disc formation could meet 
with sonic difficulties: and, ii), that entropy is conserved in reversible adiabatic processes, because 
the violation of this physical principle could produce spurious effects at galaxy scales. By the 
moment, star formation (SF) processes have to be modelled, either inspired in kpc or pc scale 
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hydrodynamical simulations (Padoan et al. 1999; Avila-Reese k, Vazquez-Semadeni 2001; Vazquez- 

Semadeni et al. 2000; Wada & Norman 2001; Kritsuk & Norman 2002) or other considerations 
(Katz 1992; Tissera, Lambas & Abadi 1997; Kennicutt 1998; Yepes el al. 1997; Silk 2001; Springel 
& Hernquist 2002b; Elmegreen 2003; Padoan et al. 2003). 

The first choice to be made when designing this kind of codes is the gravity solver. Among 
current numerical methods, those that employ adaptive techniques in regions of high density, either 
from a Lagrangian (as AP3M, Couchman 1991; Couchman, Thomas, k, Pearce 1995) or Eulerian 
(as the ART and MLAPM codes, Kravtsov, Klypin & Khokhlov 1997; Knebe, Green k Binney 
2001) approach, are the most suitable to meet the requirements of resolution and large dynamical 
range, accuracy and rapidity. A detailed comparison between AP3M and ART codes has been 
carried out by Knebe et al. (2000). They have found out that these codes produce results that 
arc consistent within a 10%, provided that A^steps/-D-R > 2 (A^gteps is the number of integration 
steps; DR stands for the dynamical range). The choice of the gravity solver in a cosmological 
simulation depends on its purpose. To study galaxy formation and evolution, Lagrangian codes 
have the advantage over Eulerian codes that they permit to go backwards and forwards in time 
in a very easy way. For example, the constituent particles of a given object can be identified at a 
given redshift, zi, and one can then analyze their positions in phase space and the properties of 
the objects or structures they form at a different redshift, Z2- This is a very convenient method 
to study evolutionary processes and it motivates our choice of an AP3M-based method as gravity 
solver for our simulations. 

To solve the hydro dynamical equations (and, in general, any hyperbolic system of equations 
in partial derivatives), there are also two basic different techniques: i), Eulerian methods, and, ii), 
Lagrangian methods. Eulerian methods are based on the so-called Godunov algorithm (Godunov 
1959). Their new formulations, using adaptive mesh refinements (Norman & Bryan 1998; Klein et 
al. 1998; Teyssier 2002), arc particularly well suited to combine with ART-likc codes when both 
gravitational and hydrodynamical forces are considered. For a comparison of the performances of 
a number of hydrodynamical codes of both kinds see (Kang et al. 1994; Prenk et al. 1999). 

Most of the lagrangian methods used in astrophysics are based on the SPH (Smooth Particle 
Hydrodynamics) technique (Lucy 1977; Gingold & Monaghan 1977; Monaghan 1992). Given the 
convenience of this technique to be applied to cosmological simulations, a number of authors de- 
veloped SPH codes to be used in a cosmological context. Some of them follow. Evrard (1988) first 
used SPH techniques in cosmological simulations (a P3M-SPH code). Hernquist & Katz (1989); 
Katz, Weinberg, & Hernquist (1996), as well as Navarro & White (1993), coupled a SPH code to 
the Barnes & Hut (1986) Tree algorithm. Vedel, Hellsten k Sommer-Larsen (1999) modelled their 
TREESPH code after Hernquist k Katz (1989) and Dave, Dubinski, k Hernquist (1997) introduced 
a parallel version of this code, while Steinmetz (1996) makes use of a special purpose hardware to 
compute the gravitational forces by direct summation (GRAPE). Serna, Alimi, k Chieze (1996) 
coupled SPH with a PP algorithm in a code designed to be run on a Connection Machine, and 
Alimi et al. (2002) incorporated in a Tree-SPH code the so-called Vh terms (see below). GADGET 
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(Springel, Yoshida & White 2001) uses either a Tree scheme or GRAPE, with individual integration 
timesteps, and both, serial and parallel versions. Another parallel Tree-Sph code is GASOLINE 
(Borgani et al. 2002). 

As already mentioned, AP3M-based codes are particularly well suited to study galaxy forma- 
tion through self-consistent cosmological simulations. The first AP3M-SPH code was introduced 
by Couchman, Thomas, & Pearce (1995) (Hydra code, see also Pearce &; Couchman 1997). Tis- 
sera, Lambas & Abadi (1997) carried out a second implementation. In these implementations, the 
integration timestep is global (i.e., at a given time, the same for all particles). In cosmological 
simulations, however, multiple time scales appear, due to their very large dynamical ranges from 
very dense volumes to very rarefied zones. To get an accurate enough integration scheme, and, at 
the same time, to avoid that particles in denser volumes slow down the simulations, it is advanta- 
geous to introduce individual integration timesteps, i.e., at each time, different timesteps for each 
particle, depending on the density of the region it samples. This is the optimal design of the code 
to increase the mass resolution. 

Another shortcoming of conventional SPH formulations concerns the entropy violation of the 
dynamical equations, related to the space dependence of the smoothing length of SPH particles, 
h{r,t), as noted by some authors (Hernquist 1993; Nelson &; Papaloizou 1993, 1994; Serna, Alimi, 
& Chieze 1996). A rigorous formulation of SPH requires that additional terms must be included in 
the particle equations of motion which account for the variability of h, usually termed as "the V/i 
terms" . Until very recently, they were considered as having a negligible effect on the global dynamics 
of systems (Gingold &; Monaghan 1982; Evrard 1988) and, therefore, SPH codes ignored such addi- 
tional terms and focused on energy conservation. Alternatively, treatments of hydrodynamics based 
on the Lagrange equations can be formulated that are well behaved in their conservation properties 
of both, energy and entropy, as that introduced recently by Springcl & Hernquist (2002a)^. The 
effects of entropy violation in SPH codes are not completely clear and they need to be analyzed in 
much more detail, specially in simulations where galaxies are formed in a cosmological framework. 
Previous works have analyzed this question in the case of the collapse of isolated objects and have 
found that, if such correction terms are neglected, the density peaks associated to central cores or 
shock fronts are overestimated at a ~30% level (Alimi et al. 2002). 

To make up for these shortcomings when dealing with problems related to galaxy formation 
and evolution, we introduce a new code, DEVA, where gravity is solved by means of an AP3M- 
like technique, and hydrodynamics with a SPH technique, with individual integration timesteps. 
The space dependence of the SPH resolution scale, h(r,t), has been taken into account, in order 
to ensure the conservative character of the equations of motion, as long as entropy and energy is 
considered. Another important particularity of DEVA is the attention paid to angular momentum 



^This paper is a reformulation of the so-called entropy formulation of SPH, where dynamical equations for the 
entropy instead of the energy were used (Lucy 1977; Benz & Hills 1987; Hernquist 1993), and where the energy 
conservation was not guaranteed 
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conservation, a key point to enable disc formation in simulations (Domfnguez-Tenreiro, Tissera & 

Saiz 1998; Saiz ct al. 2001, and references therein). Our choice has been to put the stress into 
conservation laws rather than into saving CPU time. But saving CPU time has also been one of 
our concerns, so that the code is fast enough that cosmological self-consistent simulations can be 
run on a modest computer machine. 

The paper is organized as follows: §1 is the Introduction. In §2, the SPH method is briefly 
reviewed and we present the SPH equations when the Vh terms are considered. § 3 is devoted 
to the particularities of the DEVA code, and § 4 to test whether the code integrates correctly the 
hydrodynamical and N-body equations (standard tests). In § 5 we introduce some self-consistent 
simulations run with DEVA, compare to one standard of reference for hydrodynamical simulations 
in a cosmological framework (the Santa Barbara cluster comparison project Frenk et al. 1999), and 
analyze the effects of the V/i terms in these simulations. Finally, in § 6, we give a summary of the 
work and discuss DEVA performances. 



2. The SPH Method 

2.1. Kernel estimates 

The basic idea of the SPH method (Lucy 1977; Gingold & Monaghan 1977) lies in representing 
the fluid elements by Ng particles which act as interpolation centers to determine the local value 
of any macroscopic variable /(r). In order to smooth out local statistical fluctuations, this inter- 
polation is performed by convolving the field /(r) with a smoothing (or kernel) function W. For 
example, the smoothed estimate of the local density is 

Pi^i) = H rnjW{rij, hi, hj) , (1) 

where =| r j — |, rrij is the mass of particle j, and hk is the smoothing length for particle k, 
which specifies the size of the averaging volume. 

Ideally, the individual particle smoothing lengths /i^ must be updated so that each particle 
has a constant number of neighbors Ns- By neighbors we mean those particles j with distances 
T~kj ^ 2/ifc. Such a condition can be exactly implemented by constructing, for each particle k, a list 
of its Ns nearest neighbors. The smoothing length of k is then defined to be 

hk = ]^\rk-Ykf\ , (2) 

where r^j. is the position vector of particle fc's most distant neighbor. Since each particle has its 
own h value, it is possible to find couples of particles (j, k) such that j is a neighbor of fc, but k is 



-6- 



not a neighbor of j. In these cases, it is obvious that the reciprocity principle^ is not satisfied and, 
therefore, simulations will not conserve momentum. In order to solve this problem, it is necessary 
to symmetrize the SPH equations by using, for example, averaged kernels (Hernquist Sz Katz 1989): 



Wij = W{rij,hi, hj) = - [W{rij, hi) + Winj, hj)] . 



(3) 



A first consequence of the adopted symmetrization procedure is the specific form for the kernel 
derivatives. As a matter of fact, Eq. (3) implies that Wij is a function of three variables: rij, hi 
and hj. Consequently, its gradient ViWjk is given by: 
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(4) 



The first part of Eq. (4), which does not involve derivatives of the smoothing lengths, is the usual 
symmetrized form of ViWjk- The second part, which involves derivatives of the smoothing lengths, 
arises because of the spatial and temporal variability of h. We shall refer to terms of this type as 
'V^ terms'. Most implementations of the SPH algorithm consider only the first one and neglect 
the terms. 



2.2. Hydrodynamic equations 



The motion of particle i is determined by the momentum and energy equations: 

d-Vi p 
-df = 

dui Pi dp. 



^ = af + ar--V$, (5) 



dt Pi dt 

where $i is the local gravitational potential, af is the acceleration due to pressure forces, a^**'^ 
is the acceleration due to viscosity forces, u is the specific internal energy, P = (7 — l)pu is the 
pressure (with 7 being the constant heat ratio), and Tti is the power due to non-adiabatic heating 
or cooling processes. 

A fully consistent SPH expression for pressure forces, satisfying all conservations laws (includ- 
ing entropy conservation in reversible adiabatic problems), was obtained by Nelson & Papaloizou 
(1993, 1994): 

Ff = -Y^mj^ViPj. (7) 



^The reciprocity principle states that, if at a given time the jth particle belongs to the neighbor list of the ith 
particle, then it is mandatory that, at this same time, the ith particle belongs to the neighbor list of the jth particle 
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Using Eqs. (1) and (4) to compute ViPj, one obtains 
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dr. 
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(8) 
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where ffc = Vkkf/rkkf 



On the other hand, using Eqs. (1) and (4) to compute the dpi/dt derivative appearing in Eq. 
(6), the energy equation becomes: 



dui Pi ■^-^ 

* 3 



dW{rij,hi,hj, 



dvij 



^ij ' ^ij 



dWirij,hi)^ , dW{rij,hj) ^ 



dhi 



+ 'Hi , 



(9) 



where 



rk = 



ffcfc/ • ^kkf 



Tkkf 



(10) 



Note that Eqs. (8) and (9) have been deduced by using both spatial and time derivatives 
of the SPH density as defined by Eq. (1) with the symmetrization specified in Eq. (3), because 
compatibihty with the conservation laws requires that the SPH force and energy equations are 
evaluated in consistency with the density definition. In the case of DEVA, this requirement increases 
the CPU time per integration step. In fact, since the density pi associated to a particle i depends 
on both hi and hj, for j = 1, ...,Ns, (i.e., for its Ns nearest neighbors, see Eq. 1), the computation 
of Pi at a given integration step requires the knowledge of hj for these Ns nearest neighbors at its 
beginning ^. This can be achieved either by using the hj values predicted in the previous integration 
step or by performing, at each step, a first loop over the particles to compute their hj values and, 
once it is over, a second loop to compute their hydrodynamical properties. Since we look for a high 
accuracy rather than a high computational speed, we have adopted this latter possibility. 



As usual in SPH, to account for dissipation at shocks, the above equations must be completed 
by adding an artificial viscous pressure term, Uij. When the V/i terms are considered, Hjj is added 
only to the leading term of equations (8) and (9), that is, those not involving V/i terms (Nelson & 



^As a matter of fact, each hj value must be kept fixed all along the integration step in order to avoid violating 
the reciprocity principle 
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Papaloizou 1994): 

]hi ^,^^11.. (ID 
Pi Pi ^ Pi Pj 

where we have adopted the standard viscous pressure proposed by Monaghan &; Gingold (1983): 

Pij 

where a and /? are constant parameters of order unity, rf' is a softening parameter to prevent 
numerical divergences, Cj is the local sound speed, and 



V,;ir,;i < 



P-ij 



^ijT^ij > 



3. The DEVA code 

The simulation of a system constituted by N particles usually requires a computational effort 
which considerably varies from some regions (or particles) to other. For example, regions of high 
density and submitted to strong shocks need to be siniTilatcd with timesteps much shorter than 
the rest of the system. In the AP3M+SPH codes described in the literature, all the particles in 
the system are simultaneously advanced at each timestep. The particle needing the highest time 
resolution determines the timestep length of all the others. Consequently, some few particles can 
slow down the simulation of a system. To make a code more efficient in handling with problems 
with multiple time scales, the computational effort must be centered on those particles that require 
it, avoiding useless computations for the remaining particles. In other words, it is necessary to 
allow for different timesteps for each particle. 



3.1. AP3M with individual timesteps 

A PEC (Predict-Evaluate-Correct) scheme with individual timesteps has been developed and 
implemented on our code in the following way: 

1. We enter the step n (which corresponds to the time t") with known positions r", velocities 
v", and accelerations a", for all the N particles. Furthermore, any integration scheme with 
individual timesteps needs some information to identify, at each step, those particles needing 
a recomputation of their acceleration. This information is stored in two vectors t-'*** and t"^^*, 
where t'"** is the time at which the last update of was performed, while ft^^^* = t^P-^^ + Aij 
is the time at which a recomputation of a^ will be necessary in the future. 

n n n jjiext Aast (^A\ 
i 1 ^ i 1 "-i I ^i •i''i V-*-^/ 
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2. A list is constructed with those particles j which will be advanced at the current step. Such 
particles are labelled as active. Obviously, the particle jmin with the smallest prediction time, 
j.next_fi^ must be included in this list, and fixes the timestep of the remaining active particles: 

At = mm(tf ^* - r) . (15) 

Since each step requires the update of many auxiliary arrays, it is impractical to advance only 
a single particle. For this reason, we label as active all particles within a cubic box around 
jmin- The size of the activation box is chosen, at each position, so that it contains a small 
fraction of the total number of particles. 

3. For all particles, active or not, we predict the value of r""*"-*^ and v""*"-*^ at 

f^+i = + v^At + a,"(At)V2 (16) 
v,"+i = vf + aTAi. (17) 

4. Only for active particles, we compute their accelerations and correct r"+^ and v""*"^ using : 

^„+i ^ f^+i+^(a^"+i-a^")(5t,-)V2 (18) 
v«+i = v^+i + i?(afi-a«)(<^t,), (19) 

where the choice A = 1/3 and 5 = 1/2 maintains accuracy to second order both in positions 
and velocities. In these expressions, 5tj represents the time interval elapsed from the last 
evaluation of to that performed in the current timestep 

5tj =t'' + At- tf'* . (20) 

Note that, unlike At, the 5tj value is different for each active particle. 

5. We update the global time = t" + At, as well as the t^"** and t^^^ values of each active 
particle. Here, in order to maintain the numerical stability of the AP3M algorithm, the 
individual timestep Atj must be smaller than the time scale for significant displacements or 
changes in velocity due to accelerations: 

At^ = (6Va2)V4 , (21) 
where e is the gravitational softening. 



3.2. Including SPH 

The above integration scheme may easily be extended to include hydrodynamics. The SPH 
processes involve three new independent variables in addition to those listed in Eq. (14): 

/i",<,<, (22) 
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where hj is the smoothing length, the specific internal energy, and -ii" its derivative. For all 
particles, we must then predict the value of u"'~^^ at i""*"^ 



,"+i = «" + i("At, (23) 



and compute, for active particles, both their total acceleration (a^^^^ and a^^^^^^) as well as their 
hydrodynamical variables (/i""'"^, p"""*"^, P"'+^ and ii^ ). These quantities are then used to correct 
the internal energy of active particles: 

where the choice C = 1/2 maintains accuracy to second order in internal energies. 

Now, the numerical stability requires additional limits on the timestep of each gas particle. A 
first timestep control is that concerning the time scale for significant displacements or changes in 
velocity due to accelerations: 

= {hVajfl' . (25) 

A second limit on Afj is usually given by a timestep control which combines the Courant and the 
viscous conditions: 

hi 



Atr 



_Ci + 1.2(acj + Pmaxj {/lij 



(26) 



When required, radiative cooling is implemented in an integral form (Thomas, & Couchman 
1992) using the fact that, due to the Courant condition, the density field is nearly constant over a 
time-step: 

r«.-A«r' duf"' At ^27) 



A Pi ' 

where A(T, p) is the power radiated per unit volume and Au^°°^ is the change in Ui due to cooling 
processes. 

This integral procedure circumvents the need of a control time for cooling, and, hence, it 
never limits the timestep. The numerical stability of our code requires that cooling effects must be 
updated at each step for all particles, active or not. Otherwise, the simultaneous presence of already 
cooled and not yet cooled particles in a given object would break the local pressure equilibrium 
and, as a result, cold particles would fall to the object center causing a non-physical core of very 
high density (sec §5.2 for an example). 

Fig 1 shows, for a typical cosmological simulation, the ratio of the CPU time consumed by 
an algorithm with individual timesteps to that consumed when all particles are simultaneously 
advanced. We see that the use of individual timesteps typically reduces the CPU time per step in 
a factor of five. In a pentium IV 1.7GHz personal computer, the CPU time typically consumed 
by our code in a 2 x 64^ cosmological simulation without the S/h terms is: a) 25 seconds per step 
in simulations without radiative cooling (such as the Santa Barbara cluster test of §5.1), b) from 
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25 (at high redshifts) to 70 (at low redshifts) seconds per step in cosmological simulations with 
radiative cooling (such as those of §5.2). These CPU times are increased by about 150% when the 
Vh terms are taken into account. 

4. Adiabatic Tests 

DEVA has been applied to different problems with known analytical or numerical solutions. 
The aim of such simulations was not only to test our code, but also to analyze the effects of the 
Vh terms included in it. 

4.1. The one-dimensional shock tube problem 

The one-dimensional shock tube problem proposed by Sod (1978) has become a standard test 
of all transport and source terms (including artificial viscosities) of hydrodynamic algorithms. It 
considers a perfect gas distributed on the x-axis. A diaphragm at xq initially separates two regions 
which have different densities and pressures. All particles are initially at rest. At time t = the 
diaphragm is broken and both regions start to interact. Nonlinear waves are then generated at the 
discontinuity and propagate into each region: a shock wave which moves from the high to the small 
pressure region, while the associated rarefaction wave moves in the inverse sense. The analytical 
solution to this problem has been given by Hawley, Smarr, k, Wilson (1984) and Rasio & Shapiro 
(1991). In our simulation, wc have considered N = 4096 gas particles initially distributed in the 
interval < a; < 1 according to: 

p=lP=li; = (fora:< 0.5) 

p = 0.25 P = 0.1795 v = (for x > 0.5) 

Dissipational effects, other than those associated with the artificial viscosity (with a = P = 1 and 
rj"^ = 0.01), were ignored, as well as gravitational interactions. Fig. 2 shows our results at i = 0.15. 
We see from this figure that our results are in excellent agreement with the analytical solutions. 
The resulting profiles both in the shock wave (located at x ~ 0.75) and in the contact discontinuity 
(located at x ~ 0.6) are much less rounded than in previous SPH computations (see, e.g. Monaghan 
&: Gingold 1983; Hernquist & Katz 1989; Rasio &: Shapiro 1991) as a result of having used a larger 
number of particles and, hence, a better resolution. We also note the almost complete suppression of 
post-shock oscillations in our results. These oscillations can be seen in the previous SPH simulations 
of this problem, especially in the velocity field, while no high-frequency vibrations are perceptible in 
our results. The weak blip observed in the pressure profile at the contact discontinuity {x ~ 0.6) is 
normal in SPH codes. Such non-physical blip has been explained by Monaghan & Gingold (1983) 
as due to the fact that the smoothed estimate of pressure is computed by using discontinuous 
quantities. It is then inevitable that P has some slight perturbation at the contact discontinuity, 
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but it has a negligible effect on the motion. In this test, simulations including the V/i terms gave 
exactly the same results as those neglecting such terms. 



4.2. Adiabatic collapse of a non-rotating gas sphere 

A 3D-problem usually considered to test hydrodynamical codes is that concerning the adiabatic 
collapse of a non-rotating gas sphere. This problem has been studied from a finite-difference method 
by Thomas (1987), and from SPH simulations by Evrard (1988) and Hernquist & Katz (1989). In 
order to facilitate the comparison of our results to those obtained by these authors, we have taken 
their same initial conditions: a gas sphere of radius R and total mass Mr, with density profile 

All the N = 4096 gas particles are initially at rest and have the same specific internal energy 
u = 0.05GMt/R- Units were taken so that 2G = Mt = 2R = 1. Initially far from equilibrium, 
the system collapses converting most of its kinetic energy into heat. A slow expansion follows and, 
at late times, a core-halo structure develops with nearly isothermal inner regions and the outer 
regions cooling adiabatically. We show in Fig. 3 different system profiles at end of the simulation. 
The solid line represents the numerical solution obtained when the Vh terms have been included, 
the dashed line represents the numerical solution obtained when these terms have been neglected, 
and the points represent the numerical solution obtained by Hernquist &; Katz (1989). We see 
that, although the solid and dashed lines are not exactly superposed, both solutions are coincident 
within the error bars. 



4.3. Interpretation of the influence of the terms 

We can understand why the V/i terms have a negligible effect in the two standard tests reported 
in this section. The effect of the V/i terms on the thermal energy can be expressed by a time scale 
defined by 

ts = it's) , (29) 

where is the ratio of the specific thermal energy Ui of particle i to the change in Ui due to the 
V/i terms 

= ^ • (30) 

When the time elapsed, t, is longer than this time scale, the V/i terms will produce a non-negligible 
effect on the thermal energy. This condition can be expressed as 

* dt' 
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where / represents the area contained by the l/ts{t') curve between t' = and t' = t. 

When these timescales are computed for the non-adiabatic tests reported in this paper, one 
obtains / = 0.001 (in the shock tube problem), and / = 0.85 (in the collapse of a non-rotating gas 
sphere) . The effect of the terms are then expected to be small. 

4.4. Quntifying the effects of the Vh terms 

Testing the effects of the Vh terms on hydrodynamical evolution can be better worked out 
in isentropic processes. It is not easy, however, to get such kind of processes in simulations of 
gas evolution because gas easily develops shocks where dissipation must occur. A possibility is 
considering the adiabatic expansion of a gas sphere from a situation of equilibrium. Using expansion 
rather than collapse has the advantage that shell crossing decreases substantially, so that viscous 
force terms can be removed and the evolution is isentropic. 

Such simulations as that shown in §4.2 lead at late times, t > 3, to equilibrium spheres with 
density profiles as that displayed in Fig. 3. We use the adiabatic expansion of such spheres to test 
the effects of the V/i terms in DEVA. Initial conditions were generated by performing a simulation 
as that described in §4.2, using different numbers of particles. At t = 3, we switchcd-off its viscous 
pressures (by setting a = /? = 0) to ensure that the subsequent evolution conserves the total 
entropy. The self-gravity was also switched-off at t = 3. In absence of gravitational interactions, 
this system expands fast and, at t = 3.3, its central density has decreased by a factor of pa 25. The 
evolution from t = 3 to t = 3.3 must conserve both the total energy and entropy. 

Table 1 shows the results for this series of simulations. We see that, when the Vh correction 
terms are neglected, energy is conserved very accurately but there exists a considerable violation of 
the total entropy (about 5% in the time interval we have considered). In the opposite case, when 
such correction terms are taken into account, both total energy and total entropy are conserved 
very accurately (about 0.02%). These results appear to be independent on the number of particles. 
We then conclude that the Vh correction terms cure entropy violation, allowing, at the same time, 
a very good energy conservation. 

5. Self-Consistent Cosmological Simulations 

The DEVA code is particularly well suited to numerically follow, in a cosmological context, 
the assembly from the field of primordial fluctuations of collapsed objects, such as galaxy clusters 
or galaxies. To illustrate DEVA performances in this situation, we briefly analyze some results 
of self-consistent simulations. Self-consistency means that initial conditions are set at high z as a 
Montecarlo realization of the field of primordial fluctuations (i.e., perturbations, characterized by a 
spectrum, to a given cosmological model), and then they are left to evolve according with Newton's 
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laws and the hydrodynamical equations. 

5.1. The Santa Barbara cluster test 

The Santa Barbara cluster problem was proposed by Prenk et al. (1999) to compare the results 
obtained from different codes. The formation of a X-ray cluster in a Cold Dark Matter (CDM) 
universe has been simulated using most of the hydrodynamic codes available at that time, setting 
a standard of reference to test newly proposed hydrodynamic codes. 

The initial conditions of this test correspond to a 3a peak of the density field smoothed with 
a Gaussian filter of radius ro = 10 Mpc according to the algorithm of Hoffman & Ribak (1991(@). 
The perturbation was centered on a periodic cubic region of side L = 64 Mpc. The cosmological 
scenario is a flat CDM universe with Hq = 50 km s~^ Mpc~^ for the Hubble constant; erg = 0.9 
for the present-day linear rms mass fluctuation in spherical top hat spheres of radius 16 Mpc; and 
rib = 0.1 for the baryon density (in units of the critical density). 64^ dark matter and 64^ baryon 
particles have been used with a softening length of 20 kpc. 

To test the influence of the V/i terms, two different simulations were run with DEVA. In one 
of these simulations, the V/i terms have been considered (SBGH) while in the other they have been 
neglected (SBnoGH) . In Figure 4, the density, temperature and entropy profiles of the cluster are 
plot. The stars represent the results obtained (Prenk et al. 1999) by Jenkins from a high-resolution 
SPH simulation using a parallel version of the Hydra code (Pearce Sz Couchman 1997), while the 
circles represent the results obtained by Bryan & Norman from a high-resolution adaptive mesh 
refinement shock-capturing code, SAMR, (Bryan et al. 1995; Bryan & Norman 1995). As previously 
remarked by Frenk et al. (1999), we see that the SPH and mesh results differ at the central region. 
This figure also shows the results obtained from our code both when the Vh terms are included 
(solid line) and neglected (dashed line). Error bars correspond to the standard deviation of the 
individual SPH data. 

We see that, now, our results differ slightly depending on whether the V/i terms have been 
included or not (the slope of the density profile flattens more rapidly in SBGH than in SBnoGH; the 
temperature profile is flat in SBGH and decreases within 100 kpc in SBnoGH; the entropy proflle is 
almost flat within 100 kpc in SBGH and decreases continuously in SBnoGH). Moreover, when the 
V/i terms are neglected, we obtain results that are similar to those of previous SPH simulations, 
and very close to Jenkins' results, obtained with a much higher resolution. When these terms 
are taken into account, the results are intermediate between previous SPH and grid results. This 
suggests that, at least in part, the difference between the SPH and grid results could be due to the 
non-physical entropy introduced by SPH codes. This non-physical entropy is negative (Alimi et al. 
2002) and, therefore, it produces objects with a smaller central temperature and a higher central 
density. 

Particular attention deserves the comparison of DEVA results with those obtained from the 
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entropy conserving SPH-Tree formulation by Springel & Hernquist (2002a), hereafter S-GADGET. 
In Figure 5 wc give a comparison of the Santa Barbara cluster entropy profiles obtained in the 
SBGH run and S-GADGET, kindly provided by Y. Ascasibar and G. Yepes (Ascasibar 2003). 
Both have been run with the same number of particles and gravitational resolution. We see that 
the agreement is very good within the error bars. So, both techniques compare very well in terms 
of entropy conservation. 

The differences found between SBGH and SBnoGH simulations can be understood on the basis 
of the timescale for the V/i terms (section 4.3). The time integral of l/tg (see Eq. 29) for this test 
is now larger than unity (/ = 1.2). 

5.2. Two Flat ACDM Cosmological Models 

We now report on several self-consistent simulations run in the context of fiat ACDM cosmo- 
logical models with the aim of studying galaxy assembly (see Table 2 for a summary). We have 
considered two different ACDM models with parameters whose values are consistent with their 
recent determinations from observations: simulations ACDMl (ACDM2) have = 0.65 (0.7), 
^^baryon = 0.07 (0.04), cTg = 1.18 (1.00) and h = 0.65 (0.70) (Lahav 2002; Netterfield et al. 2002; 
Spergel et al. 2003, and references therein). Both simulations share the same seed for the Monte- 
carlo realization of the initial fluctuation field, so that each object formed in one simulation has 
its counterpart in the other simulation. In each case, we have used 64^ DM particles and 64^ gas 
particles, in a periodic box of 10 Mpc comoving side. The gravitational softening is eg = 2.3 kpc, 
and the minimum allowed smoothing length, h = eg/2, as usual. For each cosmological model, 
two simulations have been run that are identical (they have exactly the same initial conditions and 
the same values of the cosmological parameters) except that in one case the V/i terms have been 
included (simulations ACDMl-H and ACDM2-II hereafter), while in the second case these terms 
have not been taken into account (simulations ACDMl-NoH and ACDM2-NoH hereafter, sec Table 
2). Note that in these simulations the cosmological volume is homogeneously sampled, in the sense 
that no resampling multimass technique has been used to study GLO assembly. This work can 
be considered as an extension of previous works in standard CDM models (Domfnguez-Tenreiro, 
Tissera Sz Saiz 1998; Saiz et al. 2001; Tissera 2001; Tissera et al. 2002), where cosmological 
simulations had been run with a different code based on a different numerical approach with fixed 
integration timesteps and particle masses (Tissera, Lambas &: Abadi 1997). Relative to these pre- 
vious works, DEVA opens the possibility of considering the effects of the V/i terms at scales of 
galaxies in self-consistent simulations. Moreover, the simulations we report here represent an im- 
provement of the baryonic mass resolution by factors of ~ 10 in the number of baryonic particles 
sampling a GLO of a given total mass. Also, the time resolution allowed by the multistep technique 
has been improved by a factor of ~ 30 in the denser areas of the box. 

Cooling has been implemented as described in § 3.4, where the cooling curve is that from Tucker 
(1975) and Bond et al. (1984) for an optically thin primordial mixture of H and He {X = 0.76, 
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Y = 0.24) in collisional equilibrium and in absence of any significant background radiation field. 

Concerning star formation, in this work we report on direct results (i.e., no spectroptiotometric) 
obtained with the simplest implementation of star formation in the code: through a parameteriza- 
tion, similar to those used by Katz (1992) and Tissera, Lambas & Abadi (1997), see Alimi et al. 
(2002) for details, based on the Jeans criterion for a collapsing region. Gas particles are turned 
into stars according with an inefficient Schmidt-law-like transformation rule (see Kennicutt 1998; 
Silk 2001), 

dpg ^ _dp* ^ _c*P£ .^2) 

dt dt tg ' ^ ' 

where c* is a dimensionless star-formation efficiency parameter, and tg is a characteristic time-scale 
chosen to be equal to the maximum of the local gas-dynamical time tdyn = (4vrG/9g)^-'^/^, and the 
local cooling time, tcool = Ui/ui. Equation (32) implies that the probability p that a gas particle 
forms stars in a time At is 

p = l- e-"*^*/*« . (33) 

As usual, we compute p at each time step for all eligible gas particles and draw random numbers 
to decide which particles actually form stars. 

Stellar feedback processes have not been explicitly considered, but a tuning of the efficiency 
parameters can mimic these feedback effects. Galaxy- like objects of different morphologies appear 
in the simulation: disk-like objects (DLOs), early-type-like objects (ETLOs) and irregular objects. 
DLOs contain gas in an extended disk, and most stars in a massive compact central concentration. 
In simulations with lower eg values (not reported here), stars form also in the disks, along arms. 
ETLOs are very poor in gas and their stellar component have relaxed regular ellipsoidal shapes. 
Irregulars have not defined shapes, and, in most cases, they are the product of a recent merger 
or interaction event. We note that GLOs formed in ACDM2-NoH and ACDM2-H tend to be of 
later type than their ACDMl-NoH or ACDMl-H counterparts, because of the lower values the 
parameters Obaryon and as take in ACDM2 simulations (Domfnguez-Tenreiro et al. 2003). First 
analyses of GLO formation and evolution in simulations run with the DEVA code are reported in 
Saiz et al. (2002a, b); more details will be given elsewhere (Saiz el al. 2003, in preparation). Here 
we focus on different aspects related with DEVA performances. 

One important issue related to disc formation in hydrodynamical simulations is specific angular 
momentum conservation at kpc scales (Sommer-Larsen &: Dolgov 2001; Steinmetz &: Navarro 1999; 
Dommguez-Tenreiro, Tissera &: Saiz 1998; Saiz et al. 2001). In Fig. 6, the specific total angular 
momentum at z = is represented versus mass for DLOs identified in ACDM2-NoH and ACDM2-H 
with V^ir(2.2i?(j) > 150 km {R^ is the disc scalelength, see Saiz et al. 2001). The specific total 
angular momentum is plot for dark haloes, jdm (open symbols), for the inner 83 per cent of the disc 
gas mass (i.e., the mass fraction enclosed by -Ropt = 3.2i?ci in a purely exponential disc), jg (filled 
symbols), and for the stellar component of the DLOs in our simulations, js (starred symbols). We 
see that, except for the three more massive objects, jg is of the same order as jdm; so that these gas 
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particles have collapsed conserving, on average, their angular momentum. Moreover, DLOs formed 
in our simulations are inside the box defined by observed spiral discs in this plot (Fall 1983). 
In contrast, js is much smaller than either jdm or jg, meaning that the stellar component in the 
central parts has formed out of gas that had lost an important amount of its angular momentum 
in catastrophic events or that had never acquired it. We note that this result is similar to those 
obtained in simulations run with the code by Tissera, Lambas & Abadi (1997), even if the codes 
are different, as explained above, and that the global cosmological models are also different. 

The smoothed estimate of the local gas density given by Eq. (2) and the ensuing formulation of 
SPH equations, symmetrized to ensure the reciprocity principle, has made possible jg conservation 
in an axisymmetric potential. This is a delicate crucial point in SPH codes. As stated in §2.2, 
its accurate implementation requires that the individual smoothing lengths must be completely 
updated at the beginning of each integration step, what increases the CPU time requirements. 
This complication cannot be avoided, however, when an accurate jg conservation is a key point in 
the physical processes under study. A different approach to saving CPU time is then necessary. 
Using different timesteps for each particle, useless computations for particles that do not require 
high time resolution are avoided. This saves considerable amounts of CPU time. Multistepping 
has allowed us to run these simulations, that have a considerable dynamical range, in a modest 
computing machine (a pentium IV 1.7GHz personal computer, see above). 

Let us now turn to the effects of the V/i terms at kpc scales. As stated in §4.3, their general 
effect is to correct the spurious negative entropy introduced in their absence by SPH codes, mainly 
at the central regions of collapsed objects. As a consequence, when the Vh terms are taken into 
account, dissipation by the gaseous component of a given GLO increases, so that they are more 
disordered or dynamically hotter at their central regions. Equilibrium is then attained with lower 
central baryon concentrations or densities, decreasing the amount of gas infall. But this is not the 
unique effect of V/i terms on mass distribution. It is a well known effect that dark matter is pulled in 
by baryons as they lose their energy and fall onto these central volumes of collapsed configurations 
(Dalcanton, Spergel & Summers 1997; Tissera Sz Dommguez-Tenreiro 1998). Decreasing the 
amount of gas infall translates, as a consequence, into a decreasing of the amount of dark matter 
at the GLO centers. As an illustration of this effect on both the gaseous and dark components, in 
Figure 7 (upper panel) we show the circular velocity curves for the most massive ETLO formed 
in ACDMl-H (thick lines) and its counterpart formed in ACDMl-NoH (thin lines). In the lower 
panel, the corresponding circular velocity curves are represented for a DLO formed in ACDM2- 
H (thick lines) and its counterpart in ACDM2-NoH (thin lines). In these Figures, r is the radial 
distance to the GLO center of mass, solid lines are the circular velocities, V^^^{r) = V^^{r) + V^^{r) , 
short-dashed lines and long-dashed lines stand for the dark matter (dm) and the baryonic (bar, 
both stellar and luminous) contributions, respectively, given by: 

with i = bar and dm. As can be clearly seen in these Figures, the central distributions of both, 
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dark matter and baryons, are different and in any case the concentrations are lower when the V/i 

terms are inchidcd. These central concentrations arc often quantitatively estimated in literature 
through the V^^^"^^ parameter (the maximum or peak circular velocity, see Courteau 1997; Saiz ct 
al. 2001). GLOs in ACDMl-H or ACDM2-H have lower V^^^""^ values than their ACDMl-NoH and 
ACDM2-NoH counterparts, respectively. This can be seen in Figure 7 for an ETLO and a DLO, but 
the behavior is general for any GLO. Other useful parameter to quantitatively characterize circular 
velocity curves is the logarithmic slope (LS), observationally defined for disc rotation curves as the 
slope of the straight line that fits Vcii{r) in log- log scale from R2.2 = 2.2i?d up to the last measured 
point in the rotation curve (Saiz ct al. 2001, and references therein; iij is the disk scalelength) . LSs 
are a measure of the GLO halo compactness at the scales of 10 - 30 kpcs. As illustrated in Figure 
7 for an ETLO and a GLO, a tendency for GLOs to be less compact also at these scales when the 
Vh terms are taken into account has been found in the simulations. 

A second effect of the V/t terms is related to the amounts of stars formed in a given GLO at 
its formation and all along its evolution until z = 0. In these simulations, many stars form at the 
shock fronts, where gas is compressed to very high densities. Softer shocks mean less star formation 
with the same efficiency parameters. To illustrate this effect, in Figure 8 the ratios of the total 
stellar masses for the 8 (7) more massive ETLOs produced in ACDMl-H (ACDM2-H) over the 
total stellar masses of their ACDMl-NoH (ACDM2-NoH) counterparts, are plot versus their total 
virial masses, Myir. We see that, except in one case, these ratios are smaller than the unity, as 
expected. 

Concerning sizes, for ETLOs we define the intrinsic or 3D cold baryon effective radius, rch,e, as 
the radius of the sphere enclosing half the total ETLO mass in cold baryons (i.e., cold gas or stars), 
-^cb,tot- This is a measure of the ETLO size at scales of the baryonic objects. In Figure 9 we plot 
rcb,c for the ACDMl-H (ACDM2-H) ETLOs in units of those of their ACDMl-NoH (ACDM2-NoH) 
counterparts. As expected, ACDMl-H (ACDM2-H) ETLOS have larger baryonic sizes than their 
ACDMl-NoH (ACDM2-NoH) counterparts, except for one case. 

Finally, shocks heat the gas, producing an extended diffuse hot component. The effect of 
the V/i terms in this case is to lower the temperature of this diffuse gaseous halo component of 
GLOs relative to the case when these terms are not considered. Histograms for the temperature 
distribution of the gaseous component of two of the most massive GLOs in ACDMl-H (thick lines) 
and in ACDMl-NoH (thin lines) are shown in Figure 10. To draw these histograms, all the gas 
particles within a sphere of radius riim, centered at the GLO center of mass, have been considered 
(riim is the radius where the curve of integrated plasma emission luminosity reaches its asymptotic 
value). As illustrated in the histograms for these two GLOs, the temperature distributions of 
the gaseous component for ACDMl-H and ACDMl-NoH objects do not differ substantially. In 
both cases, the gas is close to biphasic, but ACDMl-H objects are slightly colder within rum than 
ACDMl-NoH objects. They are also slightly more extended at rum scales, with the small excess of 
colder particles placed at the outskirts of the gaseous configurations. 
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The differences found between ACDMl-H and ACDMl-NoH or between ACDM2-H and ACDM2- 

NoH simulations can be understood on the basis of tfie timescale for the V/i terms (section 4.3). 
The time integral of 1 /ts (see Eq. 29) is now much larger than unity (/ ~ 80) , indicating that the 
V/i terms cannot be neglected in this kind of simulations. 

As remarked in §3.2, cooling implementation in multistcp codes has to be handled with caution: 
cooling processes need to be updated at each timestep for all particles. Otherwise, gas particles 
involved in shocks suffer from a spurious cooling and non-physical cores of high density appear, 
giving rise to extremely compact objects. As an example, in Figure 9 we plot the 3D cold baryon 
half-mass radii for a variant of ACDMl-NoH, termed ACDMl-MCool (see Table 2 for details), 
where only active particles at each timestep have been allowed to cool. In Figure 9 (starred 
symbols) we see that now the effective radii for the most massive objects are significantly smaller 
than in ACDMl-NoH; the difference becomes less important as the mass of the objects decreases 
and, as a consequence, the fraction of their constituent particles involved in shocks also decreases. 
Note that, as the comparison of filled points and starred symbols in Figure 9 shows, the combined 
effects of entropy violation and spurious cooling can produce unphysically very small objects, a 
factor of ten smaller, in some cases, than the values found when these effects are circumvented (i.e., 
the ACDMl-H simulation). 



6. Summary, Discussion and Conclusions 

We present DEVA, a multistep AP3M-likc-SPH code designed to study galaxy formation and 
evolution in connection with the global cosmological model, that uses a formulation of SPH equa- 
tions ensuring both energy and entropy conservation. 

Multistepping is introduced to save CPU time. In self-consistent cosmological simulations 
multiple time scales appear, due to their large dynamical ranges. To avoid that particles in denser 
zones slow down the simulation, and, at the same time, to get a properly accurate integration 
algorithm, it is then advantageous to use individual time steps. A comparison of the CPU time 

used in a self-consistent cosmological simulation when it is run with a global timestep or with a 
multistep scheme indicates that in the second case results at an equivalent level of accuracy are 
produced '^5 times faster. 

When a multistcp scheme is adopted, a delicate issue in the study of galaxy assembly in a 
cosmological context is the cooling implementation. In DEVA, as no cooling timescale is taken 
into consideration to fix the individual timestep for each particle, cooling must be calculated in a 
non-multistep way. Otherwise, particles involved in strong shocks would spuriously cool and form 
dense objects characterized by unphysically small baryon half-mass radii. 

On writting DEVA, particular attention was paid that conservation laws of physics (energy, 
entropy, momentum) are correctly implemented in the code, so that they hold at scales and under 
physical conditions relevant for galaxy assembly in a cosmological context. The usual formula- 
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tions of SPH equations focus on energy conservation and they violate, by construction, entropy 
conservation. As a consequence, a negative entropy is numerically produced in shocks that might 
(or might not) spuriously pollute the results, depending on the problem one studies. Different au- 
thors have addressed the issue of entropy violation in SPH codes (Nelson & Papaloizou 1993, 1994; 
Hernquist 1993; Springel & Hernquist 2002a) and have given different alternative formulations of 
its equations. We have implemented a formulation that considers explicitly the effects of the Vh 
terms in SPH equations (Nelson &; Papaloizou 1994). By taking advantage of the structure of the 
AP3M algorithm in the neighbor search, the implementation of V/i terms in the code is simple and 
noiseless. 

To test the relevance of entropy violation at shock locations under different physical conditions, 
we have studied problems and run simulations (namely, the adiabatic Santa Barbara cluster for- 
mation test, Frenk et al. (1999), and fully self-consistent cosmological simulations to study galaxy 
formation including cooling an star formation) using identical initial conditions and two different 
versions of DEVA, one that takes into account the Vh terms in SPH equations and one that does 
not take them into account. We show that entropy violation has consequences on the thermody- 
namical properties of the very central regions of the Santa Barbara cluster and on the structure at 
kpc scales of galaxy-like objects (GLOs) formed in simulations, but it does not have any apprecia- 
ble consequence on the standard non-cosmological tests of hydrodynamical codes. To understand 
the origins of these different behaviors, a criterion is introduced that allows to elucidate when en- 
tropy violation is expected to have appreciable consequences on the results. In standard tests, only 
moderate shocks and time scales are involved. 

As a consequence of the non-physical negative entropy numerically produced when the Vh 
terms are neglected, both GLOs and the cluster present more concentrated baryon density profiles 

(either star or gas). Concerning the Santa Barbara cluster test, when the V/i terms are not 
included, we get results that are similar to those of previous SPH simulations that focus on energy 
conservation. However, when these terms are considered, the results are intermediate between those 
SPH and grid results. For example, the temperature profile is decreasing within about 100 kpc of 
the cluster center when the simulation is run with standard SPH codes (including DEVA without 
V/i terms), it increases when it is run with grid codes and it is flat when DEVA + V/?, is used. 
We would like to note that the accuracy figure of entropy conservation obtained with DEVA + V/i 
and the entropy-conserving SPH- Tree code by Springel & Hernquist (S-GADGET, 2002a) compare 
quite satisfactorily, as indicated by the good agreement between the Santa Barbara cluster entropy 
profiles obtained with both codes. 

In cosmological simulations, negative entropy causes galaxy-like objects (GLOs) to be dynam- 
ically less hot and gas infall onto their central regions is artificially increased, causing, also, an 
increase of the amount of dark matter at the GLO centers. These results qualitatively agree with 
those obtained with S-GADGET by Springel & Hernquist (2002a) in their simulation of a ACDM 
cosmological model. No quantitative comparisons are possible by the moment, because a standard 
of comparison for self-consistent cosmological simulations is unfortunately not available. 
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An important result of this work is that the combined effects of entropy violation and multi- 
step cooling implementation in cosmological simulations can be particularly dramatic concerning 
the concentration of mass distribution in the galaxy-like objects they produce. For example, their 
baryon half-mass radii can be up to a factor of ten smaller than half-mass radii of GLOs produced 
in entropy-conserving non-multistep runs with DEVA. 

Concerning momentum conservation, we have used a formulation of SPH equations that is 
consistent with the smoothed estimate of the local gas density (Eqs (1) and (3)). Equations are 
symmetrized to ensure that the reciprocity principle holds (that is, if at a given time the jth. 
particle belongs to the neighbor list of the zth particle, then it is mandatory that, at this same 
time, the ith particle belongs to the neighbor list of the jth particle), so that momentum and 
angular momentum are conserved. The implementation of this principle in a SPH code increases 
considerably the CPU time per integration step, because a double loop on gas particles is necessary 
to evaluate smoothing lengths. To test angular momentum conservation, we have measured the 
specific angular momentum of discs formed in self-consistent simulations. It has been found that 
conservation is good enough to obtain simulated discs with observational counterparts, without any 
need of previous heating, as already Saiz et al. (2001) have shown (see also Governato et al. 2002). 

The use of a very high number of particles could ensure angular momentum conservation in 

SPH-tree codes (Governato et al. 2002). In this paper it has been shown that codes paying a 
particular attention to the implementation of conservation laws of physics at the scales of interest 
can attain a good level of accuracy in conservation laws with more limited resources. 
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TABLE 1 
Entropy and Energy Conservation 



Run 


V/? terms 




AE 


AS 


HTestl 


Excluded 


1024 


0.01% 


5.6% 


HTest2 


Included 


1024 


0.02% 


0.02% 


HTest3 


Excluded 


2048 


0.01% 


5.2% 


HTest4 


Included 


2048 


0.01% 


0.01% 


HTestB 


Excluded 


4096 


0.02% 


5.5% 


HTcstG 


luchidcd 


109G 


0.01%: 


0.02% 



TABLE 2 
Summary of ACDM Simulations 



Run 






0-8 


h 


Vh terms 


Cooling 


ACDMl-H 


0.65 


0.07 


1.18 


0.65 


Yes 


Non-multistep 


ACDMl-NoH 


0.65 


0.07 


1.18 


0.65 


No 


Non-multistep 


ACDMl-MCool 


0.65 


0.07 


1.18 


0.65 


No 


Multistep 


ACDM2-H 


0.70 


0.04 


1.00 


0.70 


Yes 


Non-multistep 


ACDM2-NoH 


0.70 


0.04 


1.00 


0.70 


No 


Non-multistep 
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Fig. 1. — Ratio of the CPU time consumed by an algorithm where all particles are simultaneously 
advanced to that consumed when individual timesteps are considered. 
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Fig. 2. — a) Thermal energy, b) entropy, c) density, and d) pressure profiles at t = 0.15 in the 
one-dimensional shock tube problem. Points represent the DEVA results, and solid lines are the 
analytical solutions 
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Fig. 3. — Final gas density, thermal energy, pressure and radial velocity profiles for the adiabatic 
collapse of a non-rotating gas sphere. Stars: Hernquist & Katz (1989) results. Solid (dashed) line: 
DEVA results when the V/i terms are (are not) taken into account, with their corresponding errors 
(see text) 
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r[Mpc] 

Fig. 4. — Density, temperature and entropy profiles in the Santa Barbara cluster test. The stars 
represent the SPH results obtained by Jenkins, while circles represent the results obtained by Bryan 
k, Norman from an adaptive mesh refinement code. The lines correspond to the results obtained 
from DEVA when the Vh terms are taken into account (solid line) and when these terms are 
neglected (dashed line) 
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Fig. 5. — Comparison of the Santa Barbara entropy profiles obtained from the entropy conserving 
versions of DEVA (sohd hne) and GADGET (symbols) 
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Fig. 6. — Specific total angular momentum versus mass for galaxy-like objects formed in ACDM2- 
NoH (triangles) and ACDM2-H (pentagons) with V2.2 > 150 km at z = 0. Open symbols: 
dark matter halos; filled symbols: the inner 83 per cent of the disc gas mass; starred symbols: 
central stellar component (see text). The solid (dashed) box encloses the region of the diagram 
corresponding to observed spirals (ellipticals) 



- 33 - 



600 




10 20 30 

r(kpc) 

300 I \ \ \ I \ \ \ \ \ 1 \ \ 1 1 




10 20 30 

r(kpc) 

Fig. 7. — Circular velocity curves for two ETLOs (upper panel) and two DLOs (lower panel). 
ETLOs have formed in ACDMl-H (thick lines) and ACDMl-NoH (thin lines); DLOs have formed 
in ACDM2-H (thick lines) and ACDM2-NoH (thin lines). Short-dashed and long-dashed lines stand 
for the dark matter and the baryonic (both stellar and luminous) contributions to the circular 
velocity (solid line) 
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Fig. 8. — The total stellar mass inside the virial radii of ETLOs formed in ACDMl-H (filled symbols) 
and ACDM2-H (open symbols) in units of their ACDMl-NoH and ACDM2-NoH counterparts 
(where the V/i terms have not been taken into account), versus their virial masses 
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Fig. 9. — The 3D cold baryon effective radii of ETLOs formed in ACDMl-H (filled circles) and 
ACDM2-H (open circles) in units of their ACDMl-NoH and ACDM2-NoH counterparts (where the 
V/i terms have not been taken into account), versus their virial masses. Starred symbols stand for 
3D cold baryon effective radii of ETLOs formed in ACDMl-MCool in units of their ACDMl-NoH 
counterparts (see Table and text) 
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Fig. 10. — Histograms for the temperature distribution of the gas particles samphng two of the most 
massive ETLOs formed in ACDMl-H (thick hnes) and for their counterparts formed in ACDMl- 
NoH (thin hnes) 



